Optimal Shape Design for Stokes Flow 
Via Minimax Differentiability* 
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Abstract. This paper is concerned with a shape sensitivity analysis of a viscous incompress- 
ible fluid driven by Stokes equations with nonhomogeneous boundary condition. The structure 
of shape gradient with respect to the shape of the variable domain for a given cost function is 
established by using the differentiability of a minimax formulation involving a Lagrangian func- 
tional combining with function space paramctrization technique or function space embedding 
technique. We apply an gradient type algorithm to our problem. Numerical examples show 
that our theory is useful for practical purpose and the proposed algorithm is feasible. 
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1 Introduction 

This paper deals with the optimal shape design for Stokes flow inside a moving domain. This 
problem is a basic tool in the design and control of many industrial devices such as aircraft 
wings, automobile shapes, boats, and so on. The control variable is the shape of the domain, 
the object is to minimize a cost function that may be given by the designer, and finally we can 
obtain the optimal shapes. 

The efficient computation of optimal shapes requires a shape calculus (see [7] ) which differs 
from its analog in vector spaces. It is necessary to make sense of shape gradient which is a basic 
tool to obtain necessary conditions and to provide us with gradient information required by the 
gradient type optimization methods. The velocity method (see J.Cea[3] and J.-P.Zolcsio[7, 18]) 
gave a precise mathematical meaning to this notion. 

Many shape optimization problems can be expressed as a minimax of some suitable La- 
grangian functional. The characterization of the change in geometric domain is obtained by 
velocity method. Finally the use of theorems on the differentiability of a saddle point (i.e., 
a minimax) of such lagrangian functional with respect to a parameter provides very powerful 
tools to obtain shape gradient by function space paramctrization or function space embedding 
(see [5]) without the usual study of the derivative of the state. 

The function space parametrization technique and function space embedding technique are 
advocated by M.C.Delfour and J.-P.Zolcsio to solving poisson equation with Dirichlct and Nue- 
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maim condition (see[7]). In our paper [8], we apply them to a Robin problem and give its 
numerical implementation. The purpose of this paper is to use lagrangian formulation and 
theorem on the differentiability of a minimax to study the shape sensitivity analysis for Stokes 
flow, and then give a gradient type algorithm with some numerical examples to prove that our 
theory could be very useful for the practical purpose. 

This paper is organized as follows. Section 2 is devoted to the statement of a shape opti- 
mization problem for Stokes flow. In section 3, we briefly recall the velocity method which is 
used for the characterization of the deformation of the shape of the domain, and we also give 
the definitions of Eulerian derivative and shape gradient. Then we include the divergence free 
condition directly into the Lagrange functional thanks to a multiplier which plays the role of 
the adjoint state associated with the primal pressure. This leads to a saddle point formula- 
tion of the shape optimization problem for Stokes equations with nonhomogeneous boundary 
condition. 

Section 4 is devoted to the computation of the shape gradient of the Lagrangian functional 
due to a minimax principle concerning the differentiability of the minimax formulation(sec[4, 5]) 
by Function Space Parametrization technique. 

In section 5, we compute the shape gradient by using such minimax principle coupling with 
Function Space Embedding technique and get the same expression obtained in section 4. 

Finally, in the last section, with the shape gradient information, we can establish a gradi- 
ent type algorithm to solve our problem, and numerical examples show the feasibility of our 
approach for different viscosity coefficients. 

Before closing this section, we introduce some notations that will be used throughout the 
paper. 

H m (D),m £ K, denotes the standard Sobolev space of order m with respect to the set 
D, where D is either the fluid domain O or its boundary T. Note that H°(D) = L 2 (D). 
Corresponding Sobolev spaces of vector- valued functions will be denoted by H m {D) N . 

Let u = (ui, U2, ■ ■ ■ , u^) an( l v = v 2i • ' ' i v d) be two vector functions of dimension d, 
and w be a scalar function. Dm denotes the Jacobian matrix of it, i.e., Dm = {djU.i)f j =l , and 
its transpose matrix is denoted by *Dit. We also have the following linear forms: 

d 

(u, v)a = J Q u ■ v dx = J a (u, v) dx = J Q Ui u< dx, tot, v e L 2 (fl) d ; 

»=x 

d 

a(Cl;u,v) = J n aDu :Dvdx = j n a djUidjVidx, Vtt, v S i? 1 (Sl) d ; 

d 

b(n-u,w)= - J n divuwdx = -J n E diUiwdx, tot G Vw £ L 2 (Q). 

i=i 

Note that the inner products in L 2 {Vl) d is denoted by (-, -)n, and the angle product (•, •} denotes 
the usual dot product of two vectors in this paper. 

2 Formulation of the problem 

Let D. be the fluid domain in ~R N (N = 2 or 3), and the boundary r = dfl be smooth 
is described by its velocity y and pressure p satisfying the Stokes equations: 

— aAy + Vp = / in CI 
divy = in Q 

y = g on r 
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where a stands for the kinematic viscosity coefficient. Let / £ H m (M. N ) N , and g 6 H m+3 / 2 (M. N ) 
(to > to be specified) be given satisfying the compatibility condition 

g-nds = 0, (2.2) 

then we know that the solution (y,p) belongs to iJ 1 (r2) JV x L 2 (fl) and even to H m+2 (Q) N x 
ff m+1 (fi) when T is of class C m + 2 by the regularity theorem (see[9, 17]). 

Our objective is to compute the first order "derivative" of the cost function 

J(Q) = l [\y({l)-y d \ 2 dx (2.3) 



with respect to the variational domain fi. The target velocity y d is fixed in i/ 1 (IR Ar ) Ar and 
given by the designer for some purposes. 



3 The velocity method and a saddle point formulation 

Domains Q don't belong to a vector space and this requires the development of shape calculus 
to make sense of a "derivative" or a "gradient" . To realize it, there are about three types of 
techniques: J.Hadamard[ll] 's normal variation method, the perturbation of the identity method 
by J.Simon[16] and the velocity method(see J.Cea[3] and J.-P.Zolesio[7, 18]). We will use the 
velocity method which contains the others. 

Let VeE' := C([0, r); V k (R N , R N )), where V k (R N , R N ) denotes the space of all fc- times 
continuous diffcrcntiable functions with compact support contained in and r is a small 
positive real number. The velocity field 

V{t)(x) = V(t, x), x e R N , t>0 
belongs to V k (M. N ,M. N ) for each t. It can generate transformations 

T t (V)X = x(t, X), t > 0, X e R N 
through the following dynamical system 



'^(t,X) = V(t,x(t)) 
x(0,X) = X 



(3.1) 



with the initial value X given. We denote the "transformed domain" T t (V)(fl) by f2t(V) at 
t > 0. 

Furthermore, for sufficiently small t > 0, the Jacobian J t is strictly positive: 

J t (x) := det|DT t (a;)| = detDT t (x) > 0, (3.2) 

where DT 4 (x) denotes the Jacobian matrix of the transformation T t evaluated at a point x € R N 
associated with the velocity field V. We will also use the following notation: DT t _1 (x) is the 
inverse of the matrix DT t (x) , *DT t _1 (a;) is the transpose of matrix DT t _1 (a;), and the Jacobian 
matrix of T t with respect to the boundary T is denoted by Wt = Jt|*DT t _1 n\. 
We now consider the solution (y t ,pt) on fi t of the problem 

-aAy t + Vp t = f in fi t 

divy t = 0, in fl t (3-3) 

y t =g on Tt (the boundary of ilt) 
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and the associated cost function 

J(fi t ) = \ [ \y t -y d \ 2 dx (3.4) 
£ Jn t 

We say that this functional has a Eulerian derivative at Q in the direction V if the limit 

lim J(^-^) 

exists. 

Furthermore, if the map 

V i ► dJ(f2; V) : E fc -> K 

is linear and continuous, we say that J is shape differ entiable at f2. In the distributional sense 
we have 

dJ(Jl;V) = (Q, V) C fc( R N jjjvyx-pfcfRNjjiV). (3.5) 

When J has a Eulerian derivative, we say that 3 is the shape gradient of J at fi. 

Now we shall describe how to build an appropriate Lagrangian functional that takes into 
account the divergence condition and the nonhomogcneous Dirichlet boundary condition. 

Given / e H l (R N ) and g € H 5 / 2 (R N ) , we introduce a Lagrange multiplier fi and a 
functional 

L(M,y,P,v,q,v) = / {aAy -Vp + f,v)dx + / divyqdx+ (y-g,fj,)ds (3.6) 
Jn Jn Jr 

for (y,p) e F(fi) x Q(Q), e x Q(fi), and /x e iT" 1 ^^ with 

y^^if 2 ^; p(n t ) = H 2 (n t ) N nH'(n t ) N ; Q(n t ) d ^H\n t ), 

and ilo = £1 as i = 0. 

Now we're interested in the following saddle point problem 

inf sup L(n,y 7 p,v,q,fi) 

(y,p)£Y(Q)xQ{0,) (t), ( j^) e p(j2)xQ(0)xif- 1 /2(r)« 

The solution is characterized by the following: 

• The state (y,p) is the solution of problem 

— aAy + Vp = / in Cl 

divy = in SI (3.7) 

y = g on T 

The adjoint state (v, q) is the solution of problem 

-aAv + Vq = Q in fi 

divv = in ft (3.8) 

v = on T; 

• The multiplier satisfies: /x = aDv n — qn, on T. 
Hence we obtain the following new functional, 

L(Q,y,p,v,q) — / (aAy — Vp + /, v) dx + / divyqdx + / (y — g, aDv n — qn) ds. 
Jn Jn Jr 
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To get rid of the boundary integral, the following identities are derived by Green formula, 

(y -9, Vvn) ds = [(y- g, Av) + D(y - g) : Bv] dx; 
Jr Jn 

/ (y - g,qn)ds= / [div(y-g)q+ (y - g,Vq)]dx. 
Jr Jn 

Thus we obtain the new Lagrangian: 

L(tt,y,P,v,q) = / (aAy - Vp + /, v) dx + a / [(y-g,Av) + D(y - g) : Dv}dx 
Jn Jn 

+ / divyqdx- / [div(y - g)q + (y - g, Vg)] da;. 
Jn Jn 



This domain integral is advantageous for the computation of shape gradient. 

Given a velocity field V 6 E 1 and transformed domain f2 t , we can easily verify 

J(^t) = , SU P G(n u y t ,p t ,v t ,qt) 

(y t ,pt)&Y{n t )xQ(n t ) (v t ,qt)£P(n t )xQ(n t ) 

where the Lagrangian is given by 

G(^t, y t ,Pt, v t ,qt) = F(n t ,y t ) + L(n t ,V»Pt, v t , q t ) 



(3.9) 



o/ \yt-yd\ 2 dx+ (aAy t - Wp t + f,v t )dx + / divy t g t dx 
z Jn t Jn t Jn t 

+a [(y t -g,Avt}+D(y t -g):~Dv t }dx 
Jn t 

- / [ div(y t - g)q t + (y t -g, Vg t )] dx. 
Jn t 



<n t 

and J (fit) was characterized by (3.4). 

The Lagrangian G(f2t, •) has a unique saddle point (y t ,p t ,v t ,qt) G Y(il t ) x Q(flt) x 
-P(fit) x Q(£lt) which is given by the following systems: 

State equations 

-aAy t + Vp t = f in fi t 
divy t = in Q t 



(3.10a) 



y t = g 



on r t 



Adjoint state equations 



-aAv t + Vg t = y t -y d 
divv t = 
v t = 



in fl t 
in fl t 
on T t ; 



Our objective is to get the limit 



dj(0) = lim 



i(*)-i(0) 



(3.10b) 



(3.11) 



where j(t) = J(Cl t ) = inf ( y t , Pt)eY( n t )xQ(a t ) sup (t , ti5t)eP(nt)xQ(nt) G(O t , y t ,pt, v t , q t ). 

Unfortunately, the Sobolev space Y(fl t ), Q(flt), and P(f2t) depend on the parameter f, so 
we need a theorem to differentiate a saddle point with respect to the parameter t, and there 
are two techniques to get rid of it: 

• Function space parametrization technique; 

• Function space embedding technique. 

In section 4 we will use the first case, and section 5 is devoted to the second case. We will find 
that both of them can derive the same expression for dJ(f2; V). 
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4 Function space parametrization 

This section is devoted to the function space parametrization, which consists in transporting the 
different quantities (such as, a cost function) defined on the variabie domain f2 t back into the 
reference domain f2 which does not depend on the perturbation parameter t. Thus we can use 
differential calculus since the functionals involved are defined in a fixed domain f2 with respect 
to the parameter t. 

We parameterize the functions in H m (Slt) d by elements of H m (Q) d through the transfor- 
mation: 

ip i-> ip o T^ 1 : H m (n) d -> H m (Cl t ) d , integer m > 0. 

where "o" denotes the composition of the two maps and d is the dimension of the function <p. 

Note that since T* and T t _1 are diffeomorphisms, it transforms the reference domain f2 
(respectively, the boundary T) into the new domain Qt (respectively, the boundary Tt of fit)- 
This parametrization can not change the value of the saddle point. We can rewrite (3.9) as 

J(fi t ) = inf sup G(r2 t ,yoT t - 1 ,poT^ 1 ,t;oT t ~ 1 , g oT t - 1 ). (4.1) 

(y,p)ey (n) xQ(si) (v,q)ep(n) x q(q) 

It amounts to introducing the new Lagrangian for (y,p, v, q) £ Y(ft) x Q(il) x P(ii) x Q(fi): 

G(t, y,p, v, q) =G(n u y o T^po T^ 1 , v o T t _1 , g o T^ 1 ). 
The expression for G(t, y,p, v, q) is given by 

G(t,y,p,v,q) = h(t) + h{t) +h{t) + I 4 (t), (4.2) 

where 

h(t) d = f \ I lyoT^-yfdx; 

h{t) d = f / (aA(y o T,- 1 ) - V(p o T^ 1 ) + /, « o X^ 1 ) dx + / div(y o Tf 1 )(q o X^ 1 ) dx; 
h(t) = a f [(yo Tf 1 - g, A(v o Tf 1 )) + D(y o X^ 1 - g) : D(« o T t -1 )] Ax; 

h(t) = - f [div(y o TjT 1 -g)(qo T^ 1 ) + (y o T^ 1 - g, V(g o T^ 1 )}] dx, 

and its saddle point is the solution of the following variational systems: 
State system (y*,P*) G Y(Sf) x Q(fi), V(V>,7t) G P(f2) x Q(f2), 

UfW o 77 1 , v ° rr 1 ) + fe(n^o 77V o t,- 1 ) = (/,•</> o r t ~V; f 3a) 

\b(n t ;y t oT t - 1 ,noT t - 1 ) = 0. 
Adjoint state system («', g 4 ) 6 x Q(f2), V(y>, r) g P(fi) x Q(fi), 

f a(f2 t ; «' o T t ~\ V ° T" 1 ) + b(n t ; <p o T t _1 , g* o T t ~ x ) - (y* o T^ 1 - y d , <p o T f _1 )n t , 
\&(f2 t ;«* oT t -\r oTf 1 ) = 0. 

(4.3b) 

By Green formula, the equivalent expression for G(t, y*,p*, v*, g*) is obtained: 
G(i,y* J p* J t>*,a*) = i / ly'oT^-y/dz 

- o(n t ; y* o T t - X , »* o Tf 1 ) - b(ilt;v* o T^.p* o Tp 1 ) + (/, »* ° Tf 1 )^ 
- b(n t ; y< o T f "\ g* o T" 1 ) + / (y f o T^ 1 - g, aB(v f o T t _1 )n - (g f o T^n) da. (4.4) 



G 



By the transformation T t , and the following two chain rule identities, 

D( V >oT t - 1 ) = (D V -[DT t ]- 1 )oT t - 1 ; 
div(<p oT- 1 ) = {Dtp : *{DT t }~ 1 ) o Tf 1 , 



we can rewrite it on n as 



G{t,y t ,p t ,v\q t ) = ^ f |y* - y d oT t \ 2 J t dx 

- / A(t)By t : Dt>* J t dx - / B(f)D«V J t dx + / (/ o T t X) J t dx 
Jn Jn Jn 

- [ B(t)Dy t q t J t dx + I (y t -goTuC(t)Dv t -^(noT t ))wtds (4.5) 
Jn Jr 

where the notation 

A(t)r : a d =' ^(D^r 1 ] : [cr(DT t ) _1 ], S(t)r d = f [r : *(DT f ) -1 ]; 
C(t)r = ariVTt)- 1 ^ o T t ), ^ = J t |*DT ( _1 n|. 

Similarly, the variational systems (4.3) become to 

State system £ Y(Q) x Q(f2), V(0,tt) e P(Q) x Q(fi), 

f /„ A(t)Dy* : J t dx + J n B(t)lW J t dx = f Q (f o T t , -0) J t dx. 
| / a B(t)Dy*7r J t dx = 0. 

Adjoint state system (t> 4 , g ( ) € P(O) x Q(O), V(<p, r) e P(fi) x Q(ft), 

f J a A(t)Dv f : Dip J t dx + j n B(t)D V g' J t dx = J a (y* - y d o T t , <p) J t dx, 
1 J n S(t)Dv*r J t dx = 0. 



(4.6a) 



(4.6b) 



Now we introduce the theorem concerning on the differentiability of a saddle point (or a mini- 
max). To begin with, some notations are given as follows. 
Define a functional 

Q : [0, t] x X x Y -> K 

with r > 0, and X, F are the two topological spaces. 
For any i £ [0,t], define 

and the sets 

X(t) = {x t eX: g(t) = su PyeY g{t 7 x*,y)} 
Y(t, x) = {y* G y : 0(f, x, y*) = sup yer 0(t, x, y)} 

Similarly, we can define dual functionals 

h(t) = sup inf G(t,x,y) 

and the corresponding sets 

Y(t) = {y t EY: h(t) = M xeX Q(t, x, y*)} 
X(t, y) = {x t eX: Q(t, x\y) = inL xeX G(t, x, y)} 
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Furthermore, we introduce the set of saddle points 

S(t) = {(x,y) £XxY: g(t) = G(t,x,y) = h(t)} 
Now we can introduce the following theorem (see [4] or page 427 of [7]): 
Theorem 4.1 Assume that the following hypothesis hold: 
(HI) S{t)^$, te[0,r]; 

(H2) The partial derivative d t G(t, x,y) exists in [0, r] for all 



0,2/) G 



U *(*) x 

te[0.r] 



u 



X(0) x |J Y(t) 
te[o,T] 



(H3) There exists a topology Tx on X such that for any sequence {t n : t n £ [0,t]} with 
lim t n = 0, there exists x° £ -^(0) and a subsequence {t nk }, and for each k > 1, there 

n/*oo 



exists x„ k £ X(t nk ) such that 

ft) 
(ii) 



(i) lim the Tx topology, 

n/*oo 



]imMd t g(t,x nk ,y) > d t g(0,x°,y), VyeY(O); 



(H4) There exists a topology Ty onY such that for any sequence {t n : t n £ [0, r]} wi£/i lim t„ = 

0, f/iere exists y° £ Y(0) and a subsequence {t nk }, and for each k > 1, i/iere exists 
y nk £ Y(t rik ) such that 

(i) lim y nfc = in the Ty topology, 

n/*oc 



(ii) 



\imsupd t G(t,x,y nk ) <d t g{0,x,y°), Vx G X(0). 



t\0 
k/ca 



Then there exists (x°,y°) £ X(0) x F(0) such that 

M0) = lim ^-gW 

= inf sup 0^(0,*, y) = 5*3(0, z°,y ) = sup inf $0(0, x,y) (4.7) 

T/izs means that (x°,y°) £ X(0) x 3^(0) is a saddle point of dtG(0, x, y). 

In order to apply Theorem 4.1 to our problem, we should verify the four assumptions (Hl)- 
(H4) below. 

First of all, Let's check (HI). Assume that the velocity field V £ E 1 . Choose r > small 
enough, such that there exists two constants aoj/?o(0 < «o < Po), 

a < Jt(= \Jt\)<Po, VtG [0,r]. 

Now we can follow the standard proof of the existence and uniqueness of solutions of Stokes 
equations (see [17]) to obtain that there exists a unique solution (y f ,p f , v f , q*) £ Y(fl) x Q(f2) x 
P(£l) x Q(f2) (p t and q l are unique up to a constant) and 

Vie[0,r], *(t) = {*',?'} 0, r(t) = {«*,«*} ^ 0. 
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Thus (HI) is satisfied. 

The next step is to verify (H2). The partial derivative of G(y*,p 4 , v*, q r ) with respect to the 
parameter t is characterized by 



\\V 1 - V d ° T t \ 2 divV t + My* - y d o T t , -By d V t 



d t G{y t , P t ,v t ,q t ) = / 

Jn 

- [ [A'tyDy* : Dw*Ji + A(t)Dy* : D«* divV t ] dx 
Jn 

- f [B'(t)Bv t p t J t +B{t)Dv t p t divV t }dx+ f [<D/V t ,«*>Jt + </oT t ,»*>divV t ]d« 
Jn Jn 

- [ [S'(t)DyY J t + S(i)DyY div V t ] dx + / <-DgV t , C(i)Dw* - g'(n o T t )> u*da: 
+ J (^-9° T t , [C(t)Dv* -q\no T t )] d\v T V t + [C'(t)Dt>* - j'DnVjiUt) w t ds. (4.8) 

where 

T(t) = (DT*)- 1 , T'(t) = -T(t)T>V t o T t ; 
B'(t)T = r : *T'(t)), C'(t)r = a.TT'{t) (noT t ); 
,A'(f)T : <r = a[(rT'(t)) : (<rT(t)) + (rT(f)) : (<rT'(t))] 
divrVi = divVf — DV t n ■ n. 

Since V 6 E 1 , t h-» V t and i i— ► DF; are continuous, we know that for all {y t ,p t , v*, q l ) G 
Y(Q) x Q(f2) x P(f2) x Q(fl), d t G(y t ,p t ,v t ,q t ) is well defined and exists everywhere in [0, r] 
provided that f,y d e if 1 (K JV ) iV and g e ff 5 / 2 (R Ar ) JV . 

To check (H3)(i) and (H4)(i), firstly we can readily show that there exists a positive constant 
c such that 

lly'llffi^r + llp*IU a (n) ^ c II/IIl2(r«)« 

||w*||ffi(n)^ + lk*IU 2 (o) < c||y* - y d \\mn)* 
Hence there exists subsequences (y tn ,p tn ), (v tn ,q tn ) and a priori (zi,si), (22,52) such that 

y*" zi, v tn z 2 , weakly in H 1 ^) 1 *; 
pt n _v Sl) gtn _s, S2 ^ weakly in L 2 (f2). 

Passing to the limit, (zx, Si) is characterized by 

f a(0; «!, V») + b(Q; -0, a x ) = (/, t/>)h, V-0 e P(fi); 
|&(fi;z 1)7 r) =0, VttgQ^), 

and (^27^2) satisfies: 

j a(0; z 2 , if) + 6(0; y>, s 2 ) = (zi - y d , <£>)n, V</?e P(fi); 
\6(fi;* a ,r) = 0, Vr e Q(O), 

By uniqueness, we obtain (z±,s\) = (y,p) and (22,^2) = , (?) , where (y,p) and (v,g) is the 
solution of (4.3a) and (4.3b) at t — 0, respectively, i.e., 

fa(fi;y,^) + &(fi;^,p) = (/,^)n, V^eP(JJ); 
1 6(fi;y,7r)=0, Vtt e Q(fi), 
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(4.10) 



and 

f a(Sl; v, <p) + 6(fi; <p, q) = (y- y d , tp)a, V^e P(fi); 
|&(fi;«,r) =0, VreQ(Sl). 

Furthermore, we can deduce the iT 1 (SI) ^ x L 2 (Vl)— strong convergence: (y tn ,p tn ) — > (y,p) and 
(ut", _> ( U;9 ), Hence (H3)(i) and (H4)(i) are satisfied for the H 2 (Q) N x ff^fi) -strong 
topology by the classical regularity theorem(see [9, 17]). Finally, assumptions (H3)(ii) and 
(H4)(ii) are readily satisfied in view of the strong continuity of (t, y,p) i— > d t G(t, y,p, v, q) and 
h-> d t G(t,y,p,v,q). 

Hence all the four assumptions are satisfied, and we have the Eulerian derivative: 



dJ(CL;V) 



\V - Vd? divV -(y- y d> Vy d V) 



dx- I div[(a*Du - ql){BgV)]dx 



- [ [A'(0)Dy : Bv + B'{0)Bvp - (DfV, v) + B'{0)Dyq] dx, (4.11) 
Jn 

where (y,p) and (v, q) are characterized by the variational system(4.9) and (4.10), respectively, 
and the notation 

A'(0)Dy : Dp = -a[(ByBV) : Dp + By : (BpDV)]; B'(0)t = -r : *BV. 

Expression (4.11) is a domain integral, and it is easy to find that the map 

V i ^ dJ(fi; V) : E 1 -> R 

is linear and continuous, i.e., J(f2) is shape differentiable. Then according to Hadamard- Zolesio 
structure theorem (see [7],Thm.3.6 and Cor.l, p. 348), there exists a scalar distribution W(r) £ 
V 1 {Y)' such that 

dJ(0; 7) = y W(r) < V, n > ds. 

Now we further characterize this boundary expression. Since (y,p,v,q) 6 H 3 (fl) N x i? 2 (S7) x 
H 3 (fl) N x H 2 (ft) provided that T is at less C 3 (see [17]), we can use Hadamard formula (see 
[7, 19]): 

4/ F(t,x)dx= [ f£( t ,x)dx+ [ F(t,x)(V,n t )dT t (4.12) 
dtJn t Jn t ot J dUt 

for a sufficiently smooth functional F : [0, r] x M. N — > KL So we can compute the partial 
derivative for G(t, y,p, v, q) with the expression (4.2) by using Hadamard formula, 

M </l|f)l l.=o = L'' y ~ y-f-Wi ^ + \J\y-yA 2 <v,n> ds-, 



(9/ 



= / (aA(-DyV)-'V(-Vp-V),v)dx + [ (aAy - Vp + f, -BvV) dx 
Jn Jn 

+ / [div{-ByV)q + divy(-Vq • V)]dx 
Jn 

+ J^[(aAy- Vp + f,v) + divyg](V,n)ds; 

= a / {(-D2/V,Av) + d/-g,A(-Di>V)) 
B(ByV) : Bv - B(y - g) : B{BvV)} dx+a J [(y-g, Av)+B(y-g) : Dv](V,ri) ds; 
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d_ 

dt 



{h(t)} =- [ [div(-DyV)q+ div(y-g)(-Wq-V) 

- (DyV, Vg) - (y - g, V(Vg • V))] dx - ^[div(y - g)q +(y-g 7 Vq)](V, n) ds. 



Since (y,p) and (v, q) are characterized by (4.9) and (4.10) respectively, we obtain the boundary 
expression for the shape gradient, 

dj(n ; v) = £ {h(t) + i 2 {t) + i 3 (t) + h(t)} 

ot t=0 
= J^^\y-y d \ 2 + aV(y-g) :Vv^V-nd S . (4.13) 

5 Function space embedding 

In the previous section, we have used the technique of function space parametrization in order 
to get the derivative of J (fit), i.e., 

J (fit) = hrf sup G(fl t ,y,p,v,q). (5.1) 

(y, P )£Y(n t )xQ((i t ) (tj,g) e p(n t )xQ(n t ) 

with respect to the parameter t > 0. This section is devoted to a different method based on 
function space embedding technique. It means that the state and adjoint state are defined 
on a large enough domain D (called a hold-all [7]) which contains all the transformations 
{fit : < t < t} of the reference domain fl for some small r > 0. 
For convenience, let D = R . Use the function space embedding, 

J(flt)= inf sup G(fl u y,V,V,Q). (5.2) 

where the new Lagrangian 

G(fi t ,y, v, v, c) = F(fi t ,y) + L(fi u y, v, v, Q) 

= lf \y-y d \ 2 dx+ [ {aAy -VP + f,V)dx+ [ divy Qdx 

^ Jn t Jsi t Jsi t 

+ a [ [{y-g,AV)+B(y-g):BV]dx- [ [div(y - g)Q + (y - g, VQ>] dx. (5.3) 
J Sit J Sit 

Since f,y d G H 1 (R N ) N , g G H 5 / 2 (R N ) N , and O t is sufficiently smooth, the unique solution 
(y t ,Pt,v t ,qt) of (3.10) belongs to H 3 (fl t ) N X (H 2 (fl t ) n L 2 (fl t )) X (H 3 (fl t ) N n ffjftfit)*) x 
(ff 2 (fi t ) n£g(n t )) instead of Y(fi t ) x Q(fi t ) x P(fi t ) x Q(fi t ). Therefore, the sets 

X = Y = H 3 (R N ) N x H 2 (R N ), 

and the saddle points S(t) = X(t) x Y(t) are given by 

X(t) = {(y,T)EX: y\ Qt = y t , V\n t = Pt} (5.4) 
Y(t) = {(V, Q) G Y : V|n, = » t> Gin, = «t} (5.5) 

Now we begin to verify the four assumptions of Theorem 4.1 . Firstly, we can always construct 
a linear and continuous extension(see [1]): 

II : H m (fl) d -> H m (R N ) d , d = loiN, (5.6) 
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and 

n t : H m (n t ) d H m (R N ) d . (5.7) 
Therefore we can define the extensions 

y t = u t y t , v t = n tPu and v t = n tVu Q t = n tqt , (5.8) 

of y t , pt, v t and qt. So (yt,Vt) £ X(t) and (Vt, Qt) £ Y(t), an d this shows the existence of a 
saddle point, i.e., S(t) ^ 0. Then (HI) is satisfied. 

To check (H2), we compute the partial derivative of the expression (5.3), 

d t G(si t ,y,v,v,Q)= [ [m(y,v)+w 2 (y,p,v,Q)](v,n t )ds t (5.9) 

where 

WW, V) = \\y- y d \ 2 + aD(y - g) : DV; 
m(y,V,V,Q) = {aAy -VP + f,V) + (y - g,«AV -VQ) + Qdivg. 

and n t denotes the outward unit normal to the boundary T t . 

By the previous choice of f,g and y d , and V £ ^(R^,!^), d t G(n t ,y,V, V, Q) exists 
everywhere in [0,r] for all (y, V, V, 9) £ I x 7. Hence (H2) is satisfied. 

For sufficiently smooth domains O and vector fields V £ P 1 (E Ar , R^), we have shown that 
{y t iP t ) (rcsp., (u*,g*)) converge to (y,p) (resp., in the H 2 x if 1 — strong topology as t 

goes to zero in the previous section. Hence 

y t ^y = Hy, and Vt^V = n« strongly in H 2 (R N ) N , 

and 

Vt -»• P = lip, and Q t —>-Q = Tlq strongly in ff 1 ^). 
by the following lemma. 

Lemma 5.1 (see[7]) For any integer m > 1, the velocity field V £ T) m (M N ,M. N ) and a func- 
tion $ £ H m (R N ), if 

y* -» y° m H m (n) -strong 

we have 

Y t -» r in H m (R N )- strong 

where Y t : = (n y* ) o T t ~ 1 . We afe o can s/iow i/iai i/ie above result also holds for the weak topology 
of H m (R N ). 

Furthermore, assumptions(H3)(i) and (H4)(i) are satisfied for the H 3 x H 2 — strong topology. 

Now let's check (H3)(ii) and (H4)(ii). Since (y,T,V,Q) £ X x Y and O t is sufficiently 
smooth, we can use Stokes' formula to rewrite (5.9) as: 

8tG(Ot,y, V, V, Q) = j div{ [Wi (y, V) + W 2 (y, P, V, Q)] V} da:, V (y , P) £ X, (V, Q) e Y. 
Now introduce the mapping 

(y,v,v,Q) ^ [Wi(y,v) + w 2 (y,r,v,Q)]v ■. x xy ^ h 1 {r n ) n 

which is linear and continuous. 
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Furthermore, by transformation T t , the mapping 

(t; 7>, V, Q) ^ 3 t G(fi t , y, P, V, G) = / div{[Wi(y, V) + W 3 (y, P, V, G)] V} o T t J t dx 

Jn 

from [0,r]xXxytol is continuous and (H3)(ii) and (H4)(ii) are verified. This completes 
the verification of the four assumptions of Theorem 4.1. 
Hence we obtain 

dJ(n;V)= inf sup dtG(Slt,y,V,V,Q)\t=o. (5.10) 

(y,v)ex(o) ( V,Q)ey(o) 

We also note that the expression (5.9) is a boundary integral on T t which will not depend on 
(y,V) and (V, Q) outside of fit, so the inf and the sup in (5.10) can be dropped, we then get 

dJ(fi; V) = j [Wi(y, v) + W 2 (y,p, v, q)}(V, n) ds 

However, y = g, p = and (2.2) imply W2{y,p, v , q) = on the boundary T. Finally we have 



dJ(n- : V) = 




We also find that the expression of Eulerian derivative obtained by function space embedding 
was the same as (4.13) which was obtained by the function space parametrization technique, 
but the second method is obviously quick. 

6 Gradient algorithm and numerical implementation 

In this section, we will give a gradient type algorithm and some numerical examples in two 
dimensions to prove that our previous methods (i.e. Function Space Parametrization & Function 
Space Embedding) could be very useful and efficient for the numerical implementation of shape 
problems. 

We describe a gradient type algorithm for the minimization of a cost function J (CI). As we 
have just seen, the general form of its Eulerian derivative is 

dJ(ri;V) = J vV-nds 

where v is given by a result like (4.13). Ignoring regularization, a descent direction is found by 
defining 

V=-vn (6.1) 

and then we can update the shape CI as 

Cl k = (Id + h k V)fl (6.2) 

where hk is a small descent step at fc-th iteration. 

There are also other choices for the definition of the descent direction. The method used in 
this paper is to change the scalar product with respect to which we compute a descent direction, 
for instance, iJ 1 (fi) 2 . In this case, the descent direction is the unique element d € iJ 1 (J7) 2 such 
that for every V E H^Cl) 2 , 

/ Dd : DVdx = dJ(Q;V). (6.3) 
Jn 
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The computation of d can also be interpreted as a regularization of the shape gradient, and 
the choice of H 1 ^) 2 as space of variations is more dictated by technical considerations rather 
than theoretical ones. 

The resulting algorithm can be summarized as follows: 

(1) Choose an initial shape f^o! 

(2) Compute the state system and adjoint state system, then we can evaluate the descent 
direction dk by using (6.3) with Q = fife; 

(3) Set fife+i = (Id — hkdk) Qk, where hk is a small positive real number and can be chosen 
by some rules, such as Armijo rule. 



Our numerical solutions are obtained under FreeFem+- 
havc solved the following minimization problem 



mini J * (y-y d ) 2 dx 



[13]. To illustrate the theory, we 

(6.4) 



(6.5) 



subject to 

— aAy + Vp = / in SI 

divy = in Q 

[y = on T; 

The domain Q is an annuli, and its boundary has two part: the outer boundary T ou t is a unit 
circle which is fixed; the inner boundary Ti„ which is to be optimized. We choose the target 
velocity y d = (y ld ,y 2d ) as follows: 



Vi d 



yWx 2 +y 2 - 0.2)( v / a; 2 + y 2 - 1) 



V* 2 



V 



j Via 



z{^Jx 2 + y 2 - 0.2){^x 2 + y 2 - 1) 
\J x 2 + y 2 



and the target inner boundary I\„ is a concentric circle with radius 0.2. We will solve this 
model problem with two different initial shapes: 

Case 1: A circle whose center is at origin with radius 0.4, i.e., x 2 + y 2 = OA 2 ; 
Case 2: A parabolic: x 2 /9 + y 2 /4 = 1/25. 

Now the initial mesh of the two cases are shown in Figure 6.1 and Figure 6.2. 





Figure 6.1: Initial mesh in Case 1 
with 292 nodes. 
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Figure 6.2: Initial mesh in Case 2 
with 226 nodes. 



We will use the mixed finite element method to solve the state system (4.9) and adjoint 
state system (4.10) on a triangular mesh, and the popular Pl-bubble/Pl finite element couple 
(see [10]) is chosen for the velocity- pressure couple. We run the program on a home PC. 




-0.25 -0.2 -0.15 -0.1 -0.05 0.05 0.1 0.15 0.2 



Figure 6.3: Case 1, a = 1, CPU time after 9 iterations: 37.235 second 




-0.25 -0.2 -0.15 -0.1 -0.05 0.05 0.1 0.15 0.2 



Figure 6.4: Case 1, a = 0.1, CPU time after 9 iterations: 36.125 second 
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In Case 1, Figure 6.3 — Figure 6.6 give the comparison between the target shape with iterated 
shape for the viscosity coefficient a = 1,0.1,0.01,0.001, respectively. We can find that for 
a = 1, 0.1, 0.01, we have nice reconstruction, but for a = 0.001, the result is not so satisfied in 
Figure 6.6. 




Figure 6.8: Case 2, a = 0.1, CPU time after 15 iterations: 64.172 second 
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Figure 6.9: Case 2, a = 0.01, CPU time after 15 iterations: 66.141 second 




Number of iteration 

;ure 6.10: Convergence history of the cost function in two cases with a — 0. 
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In Case 2, Figure 6.7 — Figure 6.9 represent the comparison between the target shape with 
iterated shape for the viscosity coefficient a = 1,0.1,0.01, respectively. It can be shown that 
for fixed viscosity, Case 1 has better reconstruction than Case 2, that's to say, the iteration 
process depends on the choice of the initial shape. 

Figure 6.10 shows the fast convergence of our cost function (6.4) in Case 1 and Case 2 for 
the viscosity a = 0.01. 

Finally, the numerical examples show the feasibility of the proposed iteration algorithm and 
further research is necessary on efficient implementations. 
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